{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lumping a mass matrix\n",
    "\n",
    "Explicit dynamics simulations require the usage of lumped mass matrices i.e. diagonal mass matrices for which inversion can be done explicitly.\n",
    "\n",
    "We show how to do this for P1 and P2 Lagrange elements. The former case can be done quite easily whereas the latter case requires a special integration scheme which is not present by default in FEniCS.\n",
    "\n",
    "## P1 elements\n",
    "\n",
    "### First method\n",
    "\n",
    "This problem has been already tackled in the past in old forum posts, [see for instance here](https://fenicsproject.org/qa/4284/mass-matrix-lumping/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistent mass matrix:\n",
      " [[0.083 0.042 0.042 0.   ]\n",
      " [0.042 0.167 0.083 0.042]\n",
      " [0.042 0.083 0.167 0.042]\n",
      " [0.    0.042 0.042 0.083]]\n",
      "Lumped mass matrix:\n",
      " [[0.167 0.    0.    0.   ]\n",
      " [0.    0.333 0.    0.   ]\n",
      " [0.    0.    0.333 0.   ]\n",
      " [0.    0.    0.    0.167]]\n"
     ]
    }
   ],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "mesh = UnitSquareMesh(1, 1)\n",
    "V1 = FunctionSpace(mesh, \"CG\", 1)\n",
    "v = TestFunction(V1)\n",
    "u = TrialFunction(V1)\n",
    "\n",
    "mass_form = v*u*dx\n",
    "mass_action_form = action(mass_form, Constant(1))\n",
    "\n",
    "M_consistent = assemble(mass_form)\n",
    "print(\"Consistent mass matrix:\\n\", np.array_str(M_consistent.array(), precision=3))\n",
    "\n",
    "M_lumped = assemble(mass_form)\n",
    "M_lumped.zero()\n",
    "M_lumped.set_diagonal(assemble(mass_action_form))\n",
    "print(\"Lumped mass matrix:\\n\", np.array_str(M_lumped.array(), precision=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For explicit dynamics simulation, the mass matrix can then be manipulated using the diagonal vector, to compute for instance $M^{-1}w$ where $w$ is some Function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_vect = assemble(mass_action_form)\n",
    "w = Function(V1)\n",
    "iMw = Function(V1)\n",
    "iMw.vector().set_local(w.vector().get_local()/M_vect.get_local())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second method\n",
    "\n",
    "The above trick consists in summing all rows and affecting the corresponding mass to the corresponding diagonal degree of freedom. E.g. for element $T$ with number of vertices $n$, the diagonal lumped mass $i$ is obtained as:\n",
    "\n",
    "$$\\begin{equation}\n",
    "\\overline{M}_i = \\sum_j M_{ij} = \\sum_j \\int_{T} N_i(x)N_j(x) \\text{ dx} = \\int_T N_i(x) \\text{ dx} = \\dfrac{|T|}{n}\n",
    "\\end{equation}$$\n",
    "\n",
    "for linear shape functions.\n",
    "\n",
    "We can obtain the same result by integrating the mass form with a \"vertex\" scheme (i.e. quadrature points are located at the vertices of the simplex with equal weights $\\dfrac{|T|}{n}$:\n",
    "\n",
    "$$\\begin{align}\n",
    "M_{ij} &= \\int_{T} N_i(x)N_j(x) \\text{ dx('vertex')} \\\\\n",
    "&= \\sum_{k=1}^n \\dfrac{|T|}{n}  N_i(x_k)N_j(x_k) =  \\sum_{k=1}^n \\dfrac{|T|}{n} \\delta_{ik}\\delta_{jk} = \\dfrac{|T|}{n} \\delta_{ij} \\label{diagonal-mass} \\tag{1}\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lumped mass matrix ('vertex' scheme):\n",
      " [[0.167 0.    0.    0.   ]\n",
      " [0.    0.333 0.    0.   ]\n",
      " [0.    0.    0.333 0.   ]\n",
      " [0.    0.    0.    0.167]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bleyerj/.local/lib/python3.6/site-packages/ffc/jitcompiler.py:234: QuadratureRepresentationDeprecationWarning: \n",
      "*** ===================================================== ***\n",
      "*** FFC: quadrature representation is deprecated! It will ***\n",
      "*** likely be removed in 2018.2.0 release. Use uflacs     ***\n",
      "*** representation instead.                               ***\n",
      "*** ===================================================== ***\n",
      "  issue_deprecation_warning()\n"
     ]
    }
   ],
   "source": [
    "M_lumped2 = assemble(v*u*dx(scheme=\"vertex\", metadata={\"degree\":1, \"representation\":\"quadrature\"}))\n",
    "print(\"Lumped mass matrix ('vertex' scheme):\\n\", np.array_str(M_lumped2.array(), precision=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## P2 elements\n",
    "\n",
    "The problem with the first method is that we obtain a singular matrix for P2 elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistent mass matrix:\n",
      " [[ 0.017 -0.003 -0.003 -0.011  0.     0.     0.     0.     0.   ]\n",
      " [-0.003  0.033 -0.006  0.    -0.011  0.    -0.003 -0.011  0.   ]\n",
      " [-0.003 -0.006  0.033  0.     0.    -0.011 -0.003  0.    -0.011]\n",
      " [-0.011  0.     0.     0.178  0.044  0.044 -0.011  0.044  0.044]\n",
      " [ 0.    -0.011  0.     0.044  0.089  0.044  0.     0.     0.   ]\n",
      " [ 0.     0.    -0.011  0.044  0.044  0.089  0.     0.     0.   ]\n",
      " [ 0.    -0.003 -0.003 -0.011  0.     0.     0.017  0.     0.   ]\n",
      " [ 0.    -0.011  0.     0.044  0.     0.     0.     0.089  0.044]\n",
      " [ 0.     0.    -0.011  0.044  0.     0.     0.     0.044  0.089]]\n",
      "Lumped mass matrix:\n",
      " [[0.    0.    0.    0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.333 0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.167 0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.167 0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.167 0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.    0.167]]\n"
     ]
    }
   ],
   "source": [
    "V2 = FunctionSpace(mesh, \"CG\", 2)\n",
    "v = TestFunction(V2)\n",
    "u = TrialFunction(V2)\n",
    "\n",
    "mass_form = v*u*dx\n",
    "mass_action_form = action(mass_form, Constant(1))\n",
    "\n",
    "M_consistent = assemble(mass_form)\n",
    "print(\"Consistent mass matrix:\\n\", np.array_str(M_consistent.array(), precision=3))\n",
    "\n",
    "M_lumped = assemble(mass_form)\n",
    "M_lumped.zero()\n",
    "M_lumped.set_diagonal(assemble(mass_action_form))\n",
    "print(\"Lumped mass matrix:\\n\", np.array_str(M_lumped.array(), precision=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can however adapt the second method by considering a generalized \"vertex\" quadrature scheme with quadrature points located at the element degrees of freedom points i.e. at vertices and facets mid-points for P2 elements. In this case, we still have $\\eqref{diagonal-mass}$. To do this we need to implement a custom quadrature scheme called `\"lumped\"` based on [the following this discussion](https://bitbucket.org/fenics-project/dolfin/issues/955/support-for-arbitrary-quadrature-rules)."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "This has been implemented in the file :download:`lumping_scheme.py` for triangles and tetrahedra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lumped mass matrix ('lumped' scheme):\n",
      " [[0.083 0.    0.    0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.167 0.    0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.167 0.    0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.167 0.    0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.083 0.    0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.083 0.    0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.083 0.    0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.083 0.   ]\n",
      " [0.    0.    0.    0.    0.    0.    0.    0.    0.083]]\n"
     ]
    }
   ],
   "source": [
    "import lumping_scheme\n",
    "\n",
    "M_lumped = assemble(v*u*dx(scheme=\"lumped\", degree=2))\n",
    "print(\"Lumped mass matrix ('lumped' scheme):\\n\", np.array_str(M_lumped.array(), precision=3))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Format de la Cellule Texte Brut",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
